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Abstract 

We describe an explicit algorithm to factorize an even antisymmetric 
N 2 matrix into triangular and trivial factors. The construction resembles 
the Crout algorithm for LU factorization. It allows for a straight forward 
computation of Pfaffians (including their signs) at the cost of iV 3 /3 flops. 
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Pfaffians play a role in statistical physics as well as in quantum field theories 
(QFT) related to particle physics. They are defined for antisymmetric even-sized 
quadratic matrices A with elements ay G C and i,j = 1, . . . , N = 2n. As basic 
definition we take 

Pf(^4) = T^j ^ S g n ( 7r ) a ^(l)^(2) a 7r(3)7r(4) ' * * d-n{N-l)-n(N) , (1) 

a sum over permutations in the symmetric group of iV elements reminiscent of the 
definition of determinants. In the path integral formulation of QFT one encounters 
Gaussian Grassmann integrals pQ for Majorana fermions of the form 

J D£e-*t TA * = Pf(A). (2) 

Here the components £j, i = l,...,N carry indices that stand for a compound 
labeling a Euclidean lattice site x and a Dirac spinor component. The identity 
follows from the rules of Grassmann integration with £)£ being a product over 
all differentials ordered such that the sign works out. In applications like [2 J the 
antisymmetric matrix A = C(p + m) is built from charge conjugation C and a 
(lattice) Dirac operator If). It may depend on other fields such as a scalar field 
(m — > m + tp{x)) for the Gross Neveu model pJJ, Majorana fermions and 
Pfaffians appear almost unavoidably in formulations of supersymmetric models, 
see [5] for an early attempt of a lattice simulation as well as [6] for more recent 
results. 

Up to a nontrivial sign a Pfaffian is the square root of a determinant. This 
may be shown by doubling the Majorana fermion §2§ into a Dirac fermion 

[Pf (A)} 2 = J L>£L>£'e-^ T ^ +e ' T ^ = J D^D^e^ = det(A) (3) 

with ift = (£ + i£')/y/2, ij) = (£ — i£') T / \/2 and a corresponding definition of 
Dip Dip. If we choose the Majorana representation for the Dirac matrices, then A is 
manifestly real. Then there is a pseudofermion representation for each degenerate 
fermion paii0 

[Pf(A)] 2 = [det(A 2 )] 1/2 = J Z^e-^^-V ( 4 ) 

in terms of real bosonic variables ip which can be a starting point of a hybrid 
Monte Carlo simulation as in [3]. 

For smaller systems, in particular in two dimensions and for algorithmic inves- 
tigations, it can be of interest to compute Pfaffians and determinants exactly in 

1 Note that this is not suitable to simulate a single Dirac fermion, if its Majorana components 
are mixed by non-real interactions. 
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simulations [7], [8] and for other purposes [2], even if practicable algorithms cost 
proportional to iV 3 . As the sign of the Pfafhan cannot be obtained from algorithms 
for determinants we find it of some interest to describe in this letteid an algorithm 
for the Pfaffian itself. Gaussian elimination schemes suitably adapted to antisym- 
metric matrices are described in [H] and are also mentioned or contained in the 
deeper layers of algorithms described in recent works requiring the computation of 
Pfaffians [6], [10], [11], |12J . The recursive factorization formulae worked out be- 
low are to our knowledge however not in the literature in this easily progammable 
explict form. 

For a nonsingular antisymmetric A there exists a factorization of the form 

A = PJP T (5) 

where P is lower triangular. The trivial antisymmetric matrix J with Pf(<7) = 1 
has antisymmetric 2x2 blocks around the diagonal. Its nonzero elements are 
enumerated by 

Jn-(-iy — — ( — 1)\ z = l,2,...JV. (6) 
We have learned ([5]) in [5J where it is attributed to [13]. In the form of our 
algorithm we below give a constructive proof for it. With the factorization given, 
we may change variables in the Grassmann integral ([2]) P T £ = r] with the Jacobian_| 
D£ = det(P)Dr) and, proving in passing another well-known relation, we obtain 

r N 

Pf(A) = det(P) / D V e~^ Tjr i = det(P) Pf( J) = \\p u (7) 

i=l 

The factorization (JSJ) will be constructed in a way similar to the Crout algo- 
rithm for the LU factorization [14\ of general matrices. The matrix P has more 
independent entries than A, but the factorization is rendered unique by setting for 
all odd i 

Pa = 1, p i+li = 0, i = 1,3, 5, ...,N- 1 (8) 

beside having pij = for i < j. 

The basic idea of the algorithm is to write out ()5]) in components in a special 
order by considering for i = 1, 3, 5 . . . the pairs 

i-2' 

Pji+i = Oij - ^2 (PikPjk+i - Pik+iPjk), j = i + 1, i + 2, . . . , N, (9) 

k=l 
i-2' 

{-Pi+ii+i)Pji = aj+ij - {pi + ikPjk+i -Pi+ik+iPjk), (10) 



k=l 



j = % + 2, % + 3, . . . , N. 



2 Based on the bachelor thesis by J.R., Humboldt university, 2009. 
3 Remember that the fermionic Jacobian is inverse to the bosonic one PQ. 
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Here the primed sums over k run over odd integers only and (jHJ) has been exploited. 
If we assume for a moment that we can always divide by Pi+u+i in (fTUj) then we 
may column-wise solve for the nontrivial Pj2,Pji,Pj4,Pj3, ■ ■ -Pnn- It is decisive to 
note that in this order we only encounter columns of P on the right hand sides of 
©, CEOD that have been computed before. 

Even for simple regular matrices the assumption Pi+u+i 7^ may not always be 
fulfilled in the steps with ffTOj) . Both for this reason and to improve the numerical 
precision a pivoting scheme is mandatory. To that end we introduce a permutation 
7r as in (pP) and replace A, P by re-arranged copies A', P' with matrix elements 

d'ij = <lir(i)Tr(j), Pij = Pn(i)j (11) 

and note that (0) is equivalent to A' = P'JP' T and thus 'covariant' under such 
a relabelling. We consider the above process now for A', P' with 7r initially set 
to the identity. Each time after completing (j^J) for some % we now determine the 
value j along the column where is maximal and denote it by j max . Before 

proceeding with (TTDT) we swap the entries 7r(z + 1) -H- vr(j max ). It is important 
to note that this modification does not invalidatcl any of the earlier uses of (JHJ) 
and ffT0|) . In this way we never divide by zero except when the whole column 
constructed via (Q vanishes, in which case one can show that A is singular. In 
all other cases we arrive at the factorization of the matrix A' = IL4n T where the 
matrix II implements the index permutation with ir. Then we obtain Pf(^4) from 

n 

Pi(A) = det(n~ 1 ) Pf(A') = sgn(vr) JJp^. (12) 

i=i 

Our PJP factorization for computing the Pfaffian requires approximately iV 3 /3 
flops (counting both multiplications and additions). The signum is known by 
counting the transpositions that went into building tt. 

We end by reporting on a brief test. We consider the Wilson fermion matrix 
on an L x L square lattice 



A w = c C 



2 ~ \ {( X ~ + (1 + l»¥x,x+fl} 

M=0,1 



(13) 



which is an antisymmetric matrix with iV = 2L 2 . The Dirac matrices are given 
in terms of Pauli matrices, C = ir 2 , C70 = T\, C71 = r 3 , fi is a unit vector in the 
positive fi direction, and the 5 symbols here incorporate antiperiodic boundary 



4 This is because we have generated complete columns whose entries are permuted in the same 
way on both sides of the equations. The column indices of earlier columns on the other hand are 
untouched as j ma x ^ i + 1. 
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Figure 1: Precision of Pf (x) and the factorization (+) defined in (|14|) as functions 
of the lattice size L associated with matrices of size N = 2L 2 . 



conditions in both directions thus making Aw nonsingular without a mass term. 
The Pfaffian for this matrix can be computed exactly by Fourier transformation 
[15] . To avoid a range overflow we have used this information to determine Co 
in f JT3|) such that Pf (Aw) = 1 holds up to roundoff errors. We have coded the 
procedure described above and have run it with L = 8, 9, . . . , 50 in standard 64 
bit precision. In Figure [T] we plot the deviations 

epf = | Pi(A w ) - 1|, bpjp = \\A' W - P'JP' T \\ (14) 

as functions of L where the matrix norm in epjp is given by the largest singular 
value. We note that apart form a general trend they sometimes are exceptionally 
small leading to some non-uniformity with L. 
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